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ABSTRACT 

We present a relativistic model for the stationary axisymmetric accretion flow of a 
rotating cloud of non-interacting particles falling onto a Kerr black hole. Based on 
a ballistic approximation, streamlines are described analytically in terms of timelike 
geodesies, while a simple numerical scheme is introduced for calculating the density 
field. A novel approach is presented for describing all of the possible types of orbit 
by means of a single analytic expression. This model is a useful tool for highlighting 
purely relativistic signatures in the accretion flow dynamics coming from a strong grav- 
itational field with frame-dragging. In particular, we explore the coupling due to this 
between the spin of the black hole and the angular momentum of the infalling matter. 
Moreover, we demonstrate how this analytic solution may be used for benchmarking 
general relativistic numerical hydrodynamics codes by comparing it against results of 
smoothed particle hydrodynamics simulations for a collapsar-like setup. These sim- 
ulations are performed first for a ballistic flow (with zero pressure) and then for a 
hydrodynamical one where we measure the effects of pressure gradients on the infall, 
thus exploring the extent of applicability of the ballistic approximation. 
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■ 1 INTRODUCTION 

' Matter falling down the potential well of a gravitating ob- 
, ject is the fundamental mechanism behind some of the most 
• powerful astrophy sical phenomena in the universe (see e.g. 

■ iFrank at al.ll2002( i. The more compact the central object, the 
, deeper into the well the matter can reach, and so the greater 

■ the quantity of potential energy available for extraction. 
However, for black holes (BHs), the most compact objects in 

■ the universe, purely radial infall is inefficient at converting 
kinetic energy into radiation since there is no har d surface at 
which decelerate the infalling gas (|Shapirdll97^ ). Tradition- 
ally, rotation of the accreting matter has been invoked (and 
also observed) as the means for providing, at least temporar- 
ily, centrifugal support to give time for different dissipative 
processes to take place and release part of the binding energy 
in the form of radiation. Gas rotation can lead to the forma- 
tion of a disc-like structure (Prcndcrgast & Burbidgc 1963) 
and, indeed, accretion discs around BHs are the most com- 
monly studied engines for explain ing astrophysical p henom- 
ena such as active galactic nuclei (|Genzel et al.|[201oh . X-ray 
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binar ies (|Kinelll995l ) and gamma-ray bursts (GRBs) (|Piranl 
l2004h . and they ma y be produced as a p ossible outcome 
of tidal d isruptions |Rosswog et al.l |2009| ) and binary co- 
alescence (|Lee et al.ll2010l V Comprehensive analyses of these 
systems require full-scale, numerical magnetohydrodynamic 
simulations in a curved spacetime, together with an accurate 
microphysical description of the dissipative processes, such 
as cooling and shock formation. Nevertheless, typically two 
ingredients play the leading roles in determining the overall 
accretion efficiency: the gravitational field generated by the 
central BH and rotation of the fluid. An examination of the 
combination of these last two features is the focus of the 
present study. 

Regarding the gravitational field, it has also been re- 
cognised that, in any realistic scenario, the central BH will 
posses some amount of intrinsic angular momentum, either 
because it was born with it or as a result of accretion of mat- 
ter w ith large angular momentum l|Bardeenlll970l : lBlandfordl 
Il987h . Therefore, for applications in which it is safe to neg- 
lect the surrounding mass-energy contribution to the over- 
all spacetime curvature, the exterior metric around a phys- 
ical BH will be well approximated by the Kerr solution. 
Any substantial value for the angular momentum of the BH 
will significantly affect the innermost region of an accre- 
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tion disc around the BH, exactly where one expects to find 
the highest densities, temperatures and luminosities. For in- 
stance, the inner radius of a Keplerian-like accretion disc 
around a maximally rotating BH is around six times closer 
to the central acc retor than it would be for a non-rotating 
BH (see e.g. Shapiro fc Teukolsk"^ll983l ). while the binding 
energy for the innermost stable circular orbit increases from 
~ 5.7% of the rest-mass energy for a non-rotating BH up 
to ~ 42% for a maximally rotating one. This increase in 
both the surface area of the emitting disc and the binding 
energy released can substantially boost the overall efficiency 
of the system. Moreover, the BH angular momentum may 
play a relevant role in launching and acceler ating a jet via 
e.g. t he Blandford-Znajek mechanism (Bland ford fc Znaiekl 
1 19771 ) and might also exert a torque on an accretion disc 
whic h happens to be tilted with respect to the BH rotation 
axis l|Bardeen fc Pettersonlll975l '). 

The present work constitutes a follow-up to the analytic 
acc retion model in t roduce d within a Newtonian framework 
bv iMendoza et all ll2009l) and for a Schwarzschild space- 
time by iTeieda et al.l (|2012l ) (referred to in the following as 
Paper H). Here we extend the general relativistic results of 
Paper H to a Kerr spacetime. Our aim in this series of papers 
has been to construct a toy model for the infall feeding an 
accretion disc around a BH, based solely on the two lead- 
ing ingredients determining the fluid bulk motion: gravity 
and rotation. The model is based on the assumptions of sta- 
tionarity, axisymmetry and ballistic motion, i.e. we assume 
that the fluid particles follow geodesic lines and neglect any 
deviation from their free-falling trajectories due to pressure 
gradients, magnetic flelds, self-gravity, radiative processes, 
etc. It is clear that these assumptions constitute an over- 
simplification of the real situation, but they allow us to give 
a useful analytic description of the streamlines and velocity 
field of the resulting flow. In addition, this will enable us 
to highlight the signatures of pure relativistic efi'ects in the 
accretion dynamics, due to either the strong gravitational 
field or frame dragging, that might be otherwise masked by 
a fully hydrodynamic treatment. Our analytic description of 
the streamlines is based on the extensive body of work on 
geodesic motion already existing in the literature (see e.g, 
[sharp 1979; Cha ndrasekhar .1983: . and references therein). 
Nevertheless, by utilizing some standard identities for Jac- 
obi elliptic functions, we provide here a novel approach for 
writing the solution for the radial and latitudinal motion of 
a timelike geodesic in Kerr spacetime in terms of a single 
analytic expression. 

There are several interesting astrophysical systems 
where the accretion regimes are reasonably well approxim- 
ated by the special conditi o ns of t he toy model. F or inst ance, 
Beloborodov fc lUariono^ (1200 ll ): iKumar et al] (|2008l ) and 



Zalamea fc Beloborodovl (|2009l ') used a ballistic descrip- 



tion for the infall feeding an accretion disc around a 
newly- formed BH, following the collaps e of a massiv e 
star in the s o-calle d coUapsar scenario (jWooslevI 119931 ). 
iKumar et al.l l|2008l) described the infall wi t hin N ewto- 
nian theory while Beloborodov fc Illarionovl (120011 ) and 
IZalamea fc Beloborodo^ (|2009D made general relativistic 
studies, first for a Schwarzschild BH and then for a Kerr 
one. These last two works were mainly focused on numer- 
ically solving for the structure of the disc and predicting 
a luminosity profile for the neutrino emission, while the 



infall was described approximately, considering parabolic- 
like energies for the incoming particles and restricting the 
analysis to boundary conditions with homogeneous dens- 
ity and small, uniform rotation rates. On the other hand, 
fuU-hydrodynamic numerical simulations for the c o llapsa r 
scenario were performed b y iLee fc Ramirez- Ruij (l2006l ): 
iLopez-Camara et all (|2009l ) and iTavlor et al.l (|201li ). For 
some of the simulations presented in those works, a fairly 
stationary regime is reached where the infalling matter is 
not deviating significantly from free-fall in a large fraction 
of the spatial domain of the simulation, but a shock front de- 
velops around the accretion disc which abruptly decelerates 
the incoming streamlines and marks the transition from an 
essentially ballistic regime to a hydrodynamical one (see Pa- 
per II for a comparison wi th one of the collapsar models in 
iLee fc Rami'rez-Ruij[200(^ . In those cases the present ana- 
lytic model might be a valuable tool for describing the infall- 
ing matter which feeds the accretion disc, thereby enabling 
the exploration, in a computationally efficient manner, of a 
wide range of (often uncertain) boundary conditions (e.g. 
rotation profile, accretion rates) before performing full-scale 
numerical simulations. 

In the present work we show comparisons between our 
analytic solutions and a series of 3D smoothed particle hy- 
drodynamics (SPH) simulations made using a version of the 
pu blicly availabl e code Gadget-2 (|Springel|l2005l ). modified 
by iTavlor et al.l (|201ll ) to include a simplified treatment 
of neutrino cooling and to account for approximated gen- 
eral relativistic effects related to the Kerr metric by using 
the second- order expansion pseudo-Newton ian potential de- 
veloped bv lMukhopadhvav fc Misral (|2003l ). The purpose of 
this comparison was twofold. Firstly, we wanted to show 
the utility of the toy model as a simple, practical test for 
numerical codes which include dynamical efi'ects of gen- 
eral relativity. Simulation results for pure ballistic motion 
(with the hydrodynamic forces being zeroed) can be directly 
and quantitatively compared with exact general relativistic 
results so as to test implementations of time-stepping al- 
gorithms, pseudo-Newtonian potentials, etc. Secondly, once 
the numerical features of particle motion within the simula- 
tion have been determined, one can then re-implement the 
hydrodynamic features within the code to investigate the 
effects of fluid flow behaviour, such as pressure, cooling and 
back-reaction from the accretion disc, on the particle tra- 
jectories. This second point can be viewed as an exploration 
of the validity of the ballistic approximation for describing 
the infall, showing the utility of the toy model itself for un- 
derstanding a wide variety of astrophysical scenarios in a 
relatively simple way. 

The paper is organised as follows. The general setup of 
the model is described in Section (2] while in Section [3] the 
velocity fleld of the accretion flow is given in terms of first in- 
tegrals of motion. An analytic description of the streamlines 
is given in Section [4] in terms of Jacobi elliptic functions. A 
simple numerical scheme for calculating the density field is 
presented in Section [S] Applications of the model for some 
particular boundary conditions are then described in Sec- 
tion [6] and compared against full- hydrodynamic simulations 
in Section [7] Finally, a general discussion and our conclu- 
sions are presented in Section|8] Unless otherwise stated, we 
use geometrized units for which c = G = 1. 
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2 MODEL DESCRIPTION 

In the present work we are constructing a model for the ac- 
cretion flow of a rotating cloud of non-interacting particles 
towards a Kerr BH of mass M and specific angular mo- 
mentum a. The model is based on the assumptions of sta- 
tionarity and axisymmetry. We denote the constant accre- 
tion rate by M, where the dot represents differentiation with 
respect to proper time r. In order to describe the overall ac- 
cretion flow we adopt the Boyer-Lindquist (BL) system of 
coordinates {t, r, 6, (j)). The metric line element then has 
the familiar form jMisner et al.lflQT^ 'l 



ds^ 



where 



^ dr^ + p^de^ 



4 aM r sin 



■ Atd<j) 



Esin^ 



(2.1) 



2,2 2 , 

■ r + a cos I 



[r + a 



A = - 2Mr + a, 
A sin^ 9. 



(2.2) 



Considering a set of non-interacting particles means, in 
practice, that we are making a ballistic treatment of the fluid 
flow and hence that the accretion dynamics are solely de- 
termined by the gravitational field of the BH. Under the bal- 
listic approximation, it is convenient to describe the whole 
gas cloud as a collection of equal mass test particles. We 
take as the boundary of our model a 'sphericalQ shell at 
r = ro from which test particles are continuously injected. 
The infalling particles end up being either incorporated into 
an infinitesimally thin equatorial disc or directly accreted in- 
side the BH horizon (located at r+ = M -I- KP - a^) . The 
analytic description of the infall does not include the disc 
itself where clearly a ballistic treatment is no longer valid; 
we shall just consider both disc and BH as passive sinks of 
particles and energy. We take the particle properties at ro 
to be given by specified distribution functions: 



no = n (ro, 6o) , 
ro = r (ro, 6^0) , 

e„ = e (ro, 6*0) , 

(1)0 = 4) (ro, 6^0) , 



(2.3) 
(2.4) 
(2.5) 
(2.6) 



where n is the particle number density (as measured in a 
co-moving reference frame) and f, 9, cj) are the radial, polar 
and azimuthal components of the four- velocity, respectively. 

We require the four distribution functions in Eqs. (|2.3|l - 
(|2.6|l to be differentiable and symmetric with respect to the 
equatorial plane. Additionally, in order to avoid streamlines 
intersecting before they reach the equatorial plane, two fur- 
ther conditions need to be fulfilled. First, we require that 
test particles do not have turning points in their polar or 
radial motion as they descend towards the equatorial plane. 



^ A surface r = const, in Korr spacotimc (in BL coordinates) 
defines an spheroid with cross section r'^ + a'^ sin d, rcosd^ in 
the R-z plane (see Eqs. l6.5l and l6.6Tl . 



Second, we require the mapping Oq ^ 6 to be non-singular, 
i.e. 



do 

Wo 



>0. 



(2.7) 



3 VELOCITY FIELD 

Within the ballistic approximation, streamlines of the accre- 
tion flow correspond to timelike geodesies in Kerr spacetime. 
Consider a test particle freely falling towards the central BH 
with four-velocity it'' = x'^. The stationarity and axisym- 
metry of the Kerr metric lead to the existence of four first 
integrals of the motion which will enable us to describe the 
velocity field of the accretion fiow (for a detailed derivation 
of the constants of motion, see ICartei|[i968, '). From the con- 
servation of rest mass, one has as a first conserved quantity 
the four-velocity modulus given as 



u''u^ = -l. (S.f) 

The other three constants of motion are the total specific 
energy E, the projection of the specific angular momentum 
onto the BH rotation axis £, and the Carter constant Q. 
Using BL coordinates, these quantities are given by 



E : 



2Mr\ . 2aMrsin^( 
2 

E sin^ 6 



2aMr sin^ 9 



= p'^e^ + f cot^ ( 



2 2 
£ a cos ( 



where, for convenience, we have introduced 



1. 



(3.2) 

(3.3) 
(3.4) 



(3.5) 



Once the boundary conditions in Eqs. (|2.3|l - (|2.6|) have 
been fixed, the conserved quantities in Eqs. (|3.2p - (|3.4|) are 
completely determined. Nonetheless, note that in general 
the conserved quantities will be functions of the initial polar 
angle 9o and, hence, vary from streamline to streamline. The 
value of to is calculated from Eq. (|3.ip as a function of ro, ^0 
and <po- Making use of the four integrals of motion, one gets 
the following system of equations determining the proper 
time evolution of the coordinates of the particle 



dr 



dr 

2d(j) 
P 3- = 



dr sin 
dr ' 



(3.6) 
(3.7) 
(3.8) 
(3.9) 



where 
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+2M [Q + {aE-£f] r - a^i 
Q = Q + ea^ cos^ 6 — 1^ cot^ 6, 
A = l~ aE sin^ e, 
[E{t' + a^) -al] /A. 



(3.10) 
(3.11) 
(3.12) 
(3.13) 

The signs in Eqs. (|3.6|) and p.7p are independent of each 
other and change whenever the test particle reaches a turn- 
ing point in its trajectory (though, recaU that in the present 
scenario, the radial coordinate has been required to de- 
crease monotonically as the particle approaches the equat- 
orial plane). Regarding the polar motion, since we have 
assumed mirror symmetry with respect to the equatorial 
plane, we can consider without loss of generality that the 
particle on which we are focusing is in, say, the northern 
hemisphere, i.e., Q < 9o < -k 12. For such a particle the polar 
coordinate increases from Qa to 7r/2. We then take the minus 
sign in Eq. (|3.6p and the plus sign in Eq. (I3.7|l . 

In Eqs. (|3.6|) - p.9|l . we already have expressions for the 
four components of the velocity field 



(3.14) 
(3.15) 
(3.16) 
(3.17) 



v7^ 

_A^aB sin^ d 
p2 sin^ e ' 



Eqs. (|3.14|) - (|3.17[) represent velocities with respect to 
the BL coordinate system as opposed to physical, locally 
measured ones. In ord er to get a local descr iption of the ve- 
locity field, we follow iBardeen et al] l| 19721 ) and introduce 
a set of locally non-rotating frames (LNRFs). Associated 
with each LNRF there is an orthonormal tetrad of four- 
vectors constituting a local Minkowskian coordinate set of 
basis vectors. If we denote with a bar coordinates with re- 
spect this local reference frame, the corresponding physical 
three-velocity field is given by 



7 



IP 
Ve 
IP ' 
pt 

^ sin 6* ' 



(3.18) 
(3.19) 
(3.20) 



where 7 is the Lorentz factor between the LNRF and the 
test particle passing by, and is given by 



7 = W1 + 



e 



E sin^ 9 ■ 



(3.21) 



Note that both sets of expressions for the velocity field 
(Eqs.Em - I3T7I and Eqs. lSTTSl - 13^20)1 are functions of the 
position (r, 8) as well as of the conserved quantities along 



each streamline, which have been determined by the initial 
position (ro, 80) ■ Therefore, to use them in practice we need 
to provide an explicit mapping from (ro, 80) to (r, 9). Such 
a mapping will be given in the next section in terms of an 
expression for the streamlines. 



4 STREAMLINES 

In this section we give an analytic solution for the ra- 
dial and latitudinal motion of test particles freely falling 
in Kerr spacetime. It is not within the scope of the present 
work to give an exhaustive discussion of the qualitative fea- 
tures of the trajectories. For a thorough di scussion o f time - 
like geodesies a round a rotatin g BH, see Wilk ins (1972); 
iBardc cn ( 197^); IChandrasekhaijll983i); iDyniiiikova (1989): 
!Froloy"fcNovikovl ([l99i) : lKraniotiJ(i2004^ : lFuiita fc Hikidal 
(|2009^ : [ Grossman et al. I (|2012l ). while an in-depth ana- 
lysis of the latitudinal a nd radial motion is g iven by 
Ide Felice fc Calvanil (|l972l ): iBicak fc Stuchhl3 (|l97d ). In the 
following, we build on these previous works and give a novel 
approach for expressing the radial and latitudinal solutions 
of timelike geodesic motion by means of a single analytical 
formula. 

Given the assumptions of stationarity and axisymmetry, 
all that is needed for completely describing a streamline of 
the present model is to consider the projection of an arbit- 
rary timelike geodesic onto the r-8 plane. For doing this, it 
is sufficient to consider Eqs. p.6p and (|3.7p and to combine 
them in the following way 



dr' 



d8' 



Bo 



(4.1) 



The solutions to both sides of the above equation can 

be ex pressed in terms of elliptic integrals (see e.g. 

iBvrd fc Fried man, ,1954 ). Let us introduce the following 
definitions 



$(r) 



dr' 



A9' 



(4.2) 



(4.3) 



where r„ and 8a, are as yet unspecified reference points in 
the particle trajectory. With these definitions we can rewrite 
Eq. (1331 as 



4>(r) - ^(ro) = ■f{8o) - *(e). 



(4.4) 



In the following, we give explicit expressions for $(r) and 
4.1 Radial solution 

Consider first Eq. (|4.2p . The procedure for solving this in- 
tegral is technically the same as in the Schwarzschild case 
(see Paper II), with the solution depending on the nature 
of the roots of TZ{r). The physical interpretation of these 
roots is clear: whenever they are real and greater than r+, 
they constitute turning points of the radial motion at which 
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r changes sign and the direction of integration for the ra- 
dial integ ral rever s es. In the f ollowing, we briefl y review the 
results of lWilkind (|l972h and iDvmnikoval () 19861 ) concerning 
the properties of the roots of 7l{r). 

Since TZ{r) is a fourth order polynomial, there are the 
following possibilities for the roots: all four are real; two 
are real while the other two form a complex conjugate 
pair; there are two pairs of complex conjugates. The first 
two cases include the possibility of multiplicity of the real 
roots. Although one can express the roots analytically in 
terms of the parameters of t he orbit {E, I and Q) (see e.g. 
lAbramowitz fc StegunlB-OTOl ). we do not give the final ex- 
pressions coming from this here because we do not find them 
particularly useful in practice. Instead, we assume that we 
have already found the four roots, either analytically or by 
means of a root finding algorithm, and write 7?.(r) as 

72.(7-) = e(r- — Ta){i — ri,)(r — r^){r — r^). (4.5) 
We label the roots in the following way: 

(i) If Ttir) has four real roots, and in order to satisfy 
the condition 7?.(r) > 0, there are two possibilities: either r 
is bracketed in between two non-negative consecutive roots 
of TZ{r), or r has a lower bound and is unbounded above. 
The latter case represents an open orbit with the largest 
positive root being the only turning point. 

In the first case, we call the roots bracketing r, and 
Vb (with < rt). In the second case we again take r„ as 
the lower bound for r and let be the negative root with 
the largest absolute value. In both cases the two remaining 
roots are denoted as and r<j (with |rc| < \rd\)- 

(ii) If TZ{r) has two real roots and a complex conjugate 
pair, we take r„ and (with |r„| < \rd\) to be the real roots 
while Vb and form the complex conjugate pair. 

(iii) If TZ{r) has two complex conjugate pairs of roots, 
we take = r^* and = with Re(ra) < Re(ri,). Note 
that this possibility is a special characteristic of Kerr space- 
time, since in the Schwarzschild case one of the roots was 
zero and hence at least one other root was real as well. 

With this way of l abelling the roots, we can now express 
the radial solution as (|Bvrd fc Friedmanlll954) 



2cn" 



^e{ra - rc)(rd - r^) 



{rb — ra){rd — r^) 
(r-d - rb)(r^ - r„) ' 



(4.6) 
(4.7) 



where cn~^(it, k^) is the inverse of the Jacobi elliptic func- 
tion cn{(fi, kr) with modulus k^. We refer intereste d readers 
to specialised literature on ellipt ic integrals (e.g. iHancoc^ 
I1917I : ICavlevllTgeil : lLawdenllT989l ') for a precise definition of 
Jacobi elliptic functions in terms of elliptic integrals. 

In the cases (ii) and (iii), where complex roots arise, 
intermediate steps in the calculation of Eq. (|4.6|) involve 
the use of complex quantities, although the final result 
${r) — $(ro) is always a real number. For alternative forms 



of the solution for the second and third cases, involving only 
explicitly real terms, see the Appendix. 

Finally, we note that the solution for the radial motion 
in Eq. 1)4. 6|) is formally identical to the solution given in Pa- 
per II for a Schwarzschild spacetime. 



4.2 Polar angle solution 

We next consider Eq. (|4.3p . From that equation it is clear 
that the latitudinal motion is restricted to those values of 6 
such that 0{6) 0. Just as in the radial case, the polar angle 
solution depends on the nature of the roots of the equation 
0(0) — 0. Whenever these roots belong to the natural do- 
main of 6, i.e. 9 € [0,7r], they constitute turning points of 
the polar motion at which 9 changes sign. These roots are 
straightforward to obtain after noticing that the equation 
0(6) — is equivalent to the following quadratic polyno- 
mial equation in cos^ 9 



2 4 
£ a cos 



9 + {Q + f ~ea 



= 0. 



(4.8) 



Let 9a be the turning point of the polar motion closest 
to the polar axis. For £ 0, 9a corresponds to the smallest 
positive root of Eq. (|4.8p . and in this case it is convenient to 
rewrite Q in terms of 6a as 



Q 



€■ cot^ ( 



2 2 

E a cos / 



(4.9) 



In general, the polar equation 0{9) =0 will have zero, 
two or four real roots in the interval of interest. The first case 
arises when £ — and Q > ~ea^, corresponding to a test 
particle sweeping the whole polar domain and periodically 
crossing the polar axis0 The second case corresponds to 
bounded polar motion, 9 £ [9a, n — 9a], with the test particle 
periodically crossing the equatorial plane. The third case 
corresponds to a test particle restricted to move within a 
single hemisphere (in this case the northern one) as 9 £ 
[9a, 9b], where 9b ^ tt/2 is the second turning point of the 
polar motion. 

As follows from Eqs. (|3.6|l and p.7|l . the polar and ra- 
dial motions are decoupled, and hence their turning points 
are in general independent of each other. Since we have as- 
sumed that the particle starts its journey from the northern 
hemisphere, we have 9a ^ 9o ^ tv /2. 

In terms of t he angle 9a, the solution for the polar in- 
tegral is given by l|Bvrd fc Friedmanlll954l ) 



cos 9a 



cd" 



cos 9 
cos 9a 



ke 



ke = \/— e a^/Q cos^ 



(4.10) 



(4.11) 



where cd ^{u, ke) is the inverse of the Jacobi elliptic func- 
tion cd((p, kg) with modulus kg. 

Note that *(6I„) = 0, while for 9 = tv/2 



^ Even though in this case the angles 9 = 0, n do not satisfy 
©(6) = 0, they still represent turning points since at those loca- 
tions the polar velocity changes sign discontinuously. When this 
happens, we take 9a = 0- 
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*(7r/2) 



cos 9a 



K{ke 



(4.12) 



where K{kg) is the complete elliptic integral of the first kind. 

Just as in the radial case, for some values of the para- 
meters e and Q, intermediate steps in the computation of 
Eq. (|4.10p might involve the use of complex quantities but 
the final result is always a real quantity. See the Appendix 
for alternative expressions for the polar solution involving 
just real quantities. 

For the non-rotating BH case (a — 0) one gets kg = 
from Eq. 1)4. 11|) . On the other hand, for a null value of the 
modulus, one has that cd(<^, 0) = cos((^), and therefore in 
this case Eq. (|4.10|) can be simplified as 



VQ 



cos 9 
cos 9a 



(4.13) 



which is the same expression as that found in Paper II (fol- 
lowing a different approach). 

4.3 Timelike geodesic 



Bringing together the results in Eqs. 1)4. 7[) and (|4.10p . we 
arrive at the following expression for the projection onto 
the r-9 plane of a timelike geodesic in Kerr spacetime: 



rtjrd - Vg) - rd{rb - r„)cn^ (g, k^) 
Td — Ta — {vb — r„)cn2 (^, k^) 



(4.14) 



with 



^e{ra - r^){rd - n) 



[$(ro) + *(6'o)-*(9)]. (4.15) 



Eqs. (|4.14|) and (|4.15|) constitute the analytic expression for 
the streamline of a gas element being accreted from (ro, 9o). 



5 DENSITY FIELD 

In order to get a formal expression for calculating the density 
field, we start from the continuity equation 



(n 



0, 



(5.1) 



where a semi-colon denotes covariant differentiation. Now 
we proceed to integrate the above expression over a four- 
volume element V, consisting of a streamline tube extending 
for an infinitesimal interval of coordinate time At. We take 
the spatial cross-section of such a streamline tube to be the 
collection of all the streamlines starting to fall-in from a 
differential area element dx^l^o at the initial shell and ending 
up at a second sphere with arbitrary radius r < ro. Denoting 
by dV the hypersurface delimiting the integrating volume 
and invoking the Gauss theorem, we have that 



have assumed stationarity, it is clear that the net particle 
flux through any closed spatial hypersuface at a given time 
t equals zero. Moreover, for the remaining mixed time-space 
hyperfaces of dV, the contraction u'^N^ will be, by construc- 
tion, different from zero just for a hypersurface oriented per- 
pendicularly to the radial direction. Hence, we have that 
Eq. (|5.2|) reduces to 



nu^7V^''Yl/i('-)ldtd6ld<?!) =0. (5.3) 

r 

Substituting into the previous equation that NI^^= S]^/ ^/g^ 
where g'"'" — A/p^, together with h^^^ = Ap^ sin^ 9, we arrive 
at 



" sin 9dtd9(. 



= 0. 



(5.4) 



Invoking once again the stationarity and axisymmetry con- 
ditions, it follows that dtod^o = dtdfji. Using this result in 
Eq. (15. 4p allows us to solve for n, getting 



n 



_ no Uq Po sin / d9 ^ 
u'' sin 9 \d9oy 



(5.5) 



Just as in the Schwarzschild case, analytically calculat- 
ing the partial derivative in Eq. (|5.5|) would be a rather in- 
volved process, whereas evaluating it numerically is a trivial 
task. We refer the reader to Paper II for a description of a 
numerical scheme for computing n. 

The requirement that there should be no early inter- 
sections of streamlines before the equator has been reached 
(see Ea. l2.7|l ensures that the expression for calculating the 
density in Eq. (|5.5p is well-defined. 



6 APPLICATIONS OF THE ANALYTIC 
MODEL 

We now illustrate our analytic model by applying it to an ex- 
ample scenario with boundary conditions consisting of mat- 
ter in uniform rotation on a uniform shell, i.e.: 



no = const., 
ro = const., 
4)0 ~ const., 
9o = 0. 



(6.1) 
(6.2) 
(6.3) 
(6.4) 



The condition in Eq. (|6.4p implies that, for every streamline, 

9a = 9o. 

Figure [T] shows six panels with the streamlines, velocity 
field and density contours for six different combinations of 
the flow parameters. The panels consist of spatial projec- 
tions onto the R-z plane, where R and z (together with 9) 
are the cylindrical coordinates associated with the BL ones 
and are deflned as 



(nu")^^ V=^d*s = j)nu''N^y^\h\dx = 0, (5.2) 
V ev 
where Ni^, is a unit vector normal to dV and h is the determ- 
inant of the induced metric on this hypersurface. Since we 



r^ + sin ( 



z — r cos f 



(6.5) 
(6.6) 



For specifying the set of model parameters in each case, we 
have used a, ro, VJ" and Vf, where the subscript e indicates 
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Figure 1. Streamlines, velocity field and density contours for six difi'erent combinations of the flow parameters. The values of the 
parameters used in each case are indicated above each panel. Each panel shows the spatial projection onto the R-z plane and the colour 
coding corresponds to the value of the logarithm of the particle number density Log(n /no), with the scale being indicated by the colour- 
coding bar at the left of each row. The arrows correspond to the V and components of the velocity field. The magnitude of the 
largest arrow is indicated at the bottom right of each panel. 
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Figure 2. The figure shows ro plotted against a (top panel) and 
£e (bottom panel)^ In the top panel, fixed values are taken for 
ro(= 20 M) and Ve{= 0.2) while in the bottom panel, fixed values 
are taken for ro(= 60 Af) and V^{= —0.2). Note how the BH spin 
couples with the angular momentum of the infalling matter and 
leads to a larger ro for a co-rotating flow and to a smaller To for 
a counter-rotating one. 



that the corresponding quantity is being evaluated at the 
equator of the shell. Note that, for the present boundary 
conditions, the set of parameters (a, ro, V/, V/) has a one-to- 
one correspondence with (a, ro, ro, 4>o) given by the inversion 
of the system of Eqs. (|3.18() - (|3.21[) . Also note that fixing this 
set of parameters specifies a family of models rather than a 
single one, since both length and density scales can still be 
arbitrarily and independently chosen. 

The radius of the outer edge of the disc formed as mat- 
ter reaches the equatorial plane, ro, can be calculated from 
Eq. (|4.15|l . taking first 6 — 7r/2 and then 6o = 7r/2, giving 



\/ e(ra — rc)(rd — rb) 



$(ro) 



(6.7) 



and then substituting the result into Eq. (|4.14p . 

In Figure [21 we have plotted ro, first as a function of 
the BH spin a and then as a function of the angular mo- 
mentum at the equator of the shell, i^. = l(n/2). Here we 
have assumed > 0, and so a negative value of a implies 
a counter-rotating disc. From this figure we can clearly see 
how the BH spin couples with the angular momentum of 
the disc (through the frame dragging effect), giving rise to 
a larger ro for a co-rotating disc and a smaller ro for a 



-0.8- 



a/M= 
a/M =0.5 
a/M= 1 



0.0 



0.4 



0.6 



OA 



1.0 



Figure 3. The pairs of velocity values V^-V* leading to a disc 
radius such that ro S [r+, ro] are plotted for a fixed value of 
ro = 10 A/ and the three values of a/M = 0, 0.5, 1. The upper 
boundary for each value of a represents the points {V^, V^) such 
that ro = ro, while the lower one represents those such that 
ro = r^.. Note how the domain of values in the velocity space 
leading to physically relevant accretion models shifts to smaller 
values of as a increases, because of the frame-dragging effect. 



counter-rotating one. It is also clear that, as intuitively ex- 
pected, rr) is a monotonically increasing function of both 
and VJ". 

Working with the LNRF velocities VJ and V/ makes the 
exploration of the parameter space easier since, being phys- 
ical velocities, they are naturally bounded as VJ" G (—1,0] 
and Vf G [0, 1). Furthermore, for fixed values of ro and a, a 
pair of velocities in the Ve-V* plane is also restricted by the 
condition that the resulting ro should satisfy ro G (r+,ro). 
In Figure 13] we have plotted the regions on the velocity space 
which lead to an outer radius of the disc satisfying this cri- 
terion. The plot has been constructed for a fixed value of 
ro = 10 M and three different values of a. From this fig- 
ure we observe that the domain of values in the velocity 
space leading to physically relevant accretion models shifts 
to smaller values of V/ as a increases. This behaviour is a 
consequence of the frame dragging effect: for a given test 
particle with fixed azimuthal velocity Vf, its associated an- 
gular momentum is an increasing function of a, and hence 
points in the Vg-Vf plane which, in the low-a case, did not 
have large enough angular momentum to keep the outer edge 
of the disc outside the event horizon, are able to do so for a 
larger value of a. Conversely, low-a models with an angular 
momentum only just small enough to form any disc inside 
their initial shell would have discs entirely outside their ini- 
tial shell when a is increased (thus becoming excluded from 
the parameter domain) . Also note that this parameter-space 
effect is greater on the lower boundary of Vf than on the 
upper one, which is simply due to the fact that the frame 
dragging increases as r — >■ r+. 
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7 COMPARISON WITH NUMERICAL 
SIMULATIONS 

In this section we compare the analytic solution derived 
above against SPH simulation s performed wi th the modified 
version of the code Gadget 2 (|Springe]||2od5l ) presented and 
used bv lTavlor et alj (|201ll '). In that work, the authors stud- 
ied numerically the production of the progenitors for long- 
duration GRBs as the aftermath of the collapse of a massive 
star, starting off from realistic initial conditions and includ- 
ing cooling by neutrino emission. They also included a rough 
approximation of Kerr metric general relativistic effects by 
using the second-o rder expansion pseudo-Newtonia n poten- 
tial developed by iMukhopadhvav fc Misral (|2003l ') (MM). 
Here, we compare our ballistic toy model (calculated with 
the full Kerr metric) with results from various versions of the 
SPH code implementing approximate gravity prescriptions: 
in addition to the MM potential, we have also experimented 
with the classical Newtonian po tential and the wide l y-used 
pseudo-Newtonian potential of iPaczvriskv fc Wiital (|l980l ) 
(PW). In connection with the two pseudo-Newtonian po- 
tentials, we should stress that they have been designed par- 
ticularly for capturing relevant relativistic features of im- 
portance for accretion discs, including getting correct loc- 
ations for the innermost stable circular orbit in Kerr and 
Schwarzschild spacetimes, respectively. This does not at all 
guarantee that they would be good for other purposes such 
as the infall calculations being discussed here. However, they 
have been widely used in more general contexts and so it is 
relevant for us to test them against the toy model. 

In order to perform a systematic analysis in which we 
are able to distinguish between hydrodynamic and gravita- 
tional effects, we have considered two kinds of simulation: 

(i) Ballistic free-fall, with the SPH particles being auto- 
matically removed when they cross either the equatorial 
plane or the BH horizon. Here we consider an equation of 
state (EoS) for which the fluid pressure P = 0, and hence, 
effectively 'turn off' the hydrodynamical forces. The aim of 
this kind of simulation is to highlight the differences in the 
flow dynamics coming from the use of different gravity de- 
scriptions (full general relativity, Newtonian gravity, and the 
MM and PW pseudo-Newtonian potentials) and from differ- 
ent numerical implementations of the equations of motion 
for the particles. 

(ii) FuU-hydrodynamical simulations, including back 
reaction from a growing equatorial disc and cooling in re- 
gions where the gas gets very hot (^ lO^K). Here we do not 
remove SPH particles when they reach the equatorial plane 
but rather let them settle down by themselves into a disc 
structure. In this case we employ a polytropic EoS of the 
form P = (r — 1) nil, where n is the baryon number dens- 
ity, u the internal energy per baryon and F is the polytropic 
index. We take F = 4/3, and the value of the internal en- 
ergy (taken to be constant in the initial shell) is set at an 
arbitrary but non- negligible value of one tenth of the sum 
of kinetic and potential energies of an SPH particle at the 
equator of the shell. 

For both types of simulation, we take stationary bound- 
ary conditions with SPH particles being continuously injec- 



ted with constant density and velocity distributions from a 
flxed injection radius ro- We treat the BH horizon (located 
at r+) as an inner boundary at which particles are extrac- 
ted from the simulation. For a fair comparison with the toy 
model, we report here late-time snapshots of the simulations 
in which the system has evolved to a quasi-stationary situ- 
ation (at least in the region away from the disc). As men- 
tioned above, the number of particles being used in these 
was continuously changing, but was around 2.5 x 10^ at the 
time shown. Moreover, in order to reduce the noise level 
and exploiting the axisymmetry of the system, the results 
presented in the following were obtained after averaging over 
24 cross-sectional <?!> — const, slices of the 3D simulations. 

7.1 Example I 

Here we consider a set of parameter values for the system 
that might arise in the context of collapsing stellar cores 
leading to long GRBs, namely 



a = 0.5M, (7.1) 

M = 4M0, (7.2) 

M = 0.01 Mq/s, (7.3) 

ro = 100 Af, (7.4) 

ro = -I/VSO, (7.5) 

ro0o = 0.038, (7.6) 

^0 = 0. (7.7) 



Note that, for convenience, we have used standard (non- 
geometrized) units to express the total accretion rate in 
Eg. (1731) . 

For this set of boundary conditions, in Figure |4] we 
present the analytic solution next to the results of both the 
ballistic and polytropic simulations. The figure shows a spa- 
tial projection of each case onto the R-z plane of isodensity 
contours, streamlines and velocity fields. Let us focus first on 
the ballistic simulation result (middle panel), which rapidly 
reached a stationary state. In this figure we see an overall 
satisfactory agreement with the analytic solution, although 
a closer inspection of the streamlines reveals some quantit- 
ative differences. We also note that the simulation isodensity 
contours are somewhat 'noisy' compared with the analytic 
results. Nevertheless, this level of fiuctuation is consistent 
with the effects of discretisation and interpolation within 
SPH simulations. 

FigureOshows a closer comparison of the streamlines of 
the analytic solution with the ones of the ballistic simulation. 
Since hydrodynamical effects are absent in this case, the dif- 
ferences between the numerical simulation and the toy model 
are reasonably attributed to the different descriptions of 
gravity: Kerr spacetime against the MM pseudo-Newtonian 
potential. From this figure we see that the streamlines in 
the two cases deviate importantly from each other just for 
r < 10 M and that far away from the central BH the dif- 
ferences between these two descriptions of gravity become 
negligible (bear in mind that the streamlines originate from 
r = 100 Af). 

We now investigate the inclusion of hydrodynamic prop- 
erties of the flow by comparing the polytropic SPH simula- 
tion in the bottom panel of Figure [4] with the corresponding 
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Figure 4. This figure sliows isodensity contours, streamlines and 
velocity fields for the analytic solution, and for the ballistic and 
polytropic SPH simulations for an accretion flow onto a rotating 
BH with a = 0.5 M. The model parameters are given in Eqs. 1)7. 
117.71 1. The common scale for the density colour coding is shown 
at the bottom of the figure. The velocity field in each panel is 
represented by the two-vectors (V^, V^); the length scale for 
these vectors is given at the bottom right corner of each panel. 
The SPH simulations used a varying total particle number, but 
typically this was around 2.5 X 10^ at the times shown (mass per 
SPH particle Ri 3.6 X IO-^^Mq). 




Figure 5. Streamlines corresponding to the analytic solution and 
to the ballistic SPH simulation presented in the top and middle 
panels of Figure |4] The figure shows just the first quadrant of the 
R-z plane with the BH horizon (located at r+) indicated with a 
dashed-linc quarter-circle. 




Figure 6. Streamlines corresponding to the ballistic and poly- 
tropic SPH simulations presented in the middle and bottom pan- 
els of Figure [4] The figure shows a zoom-in of the first quadrant 
of the R-z plane. The BH horizon (located at ■r+) is indicated 
with the dashed-linc quarter-circle. 



ballistic flow in the middle one. In Figure|6]we present a dir- 
ect comparison of the streamlines for these two cases. The 
first important thing to note is that in the polytropic sim- 
ulation, unlike in the ballistic one, SPH particles are not 
removed from the simulation at the equatorial plane and, 
therefore, the build-up of a disc can take place. This is in- 
deed the case for the present set of parameter values, for 
which we observe a disc that keeps growing in mass and ex- 
pands horizontally. The material in this disc corresponds to 
the fraction of the infalling matter which possesses enough 
angular momentum to remain in a stable orbit around the 
BH. Since it is not within the scope of the present article 
to study the evolution of such a disc, we show here a snap- 
shot at a time at which any kind of initial transient related 
to the initial conditions has faded away, but, at the same 
time, neither the mass nor the extension of the disc have 
grown importantly (additionally, the presence of cooling in 
the simulation aids in limiting the disc height). 

A second important feature characterising the poly- 
tropic simulation is the existence of a shock front surround- 
ing the disc and marking the boundary between two differ- 
ent flow regimes. In the pre-shock region, a clear stationary 
regime is rapidly reached where the flow moves superson- 
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PW potential 



Figure 7. Streamlines corresponding to the analytic solution and 
to three ballistic SPH simulations for a non-rotating BH (a = 0) 
while the remaining boundary conditions are as in Eqs. | |7. 211 - 117. 711 . 
The figure shows a zoom-in of the first quadrant of the R-z plane. 
The BH horizon (located at r+ and which represents the inner 
boundary) is indicated with the dashed-line quarter-circle. 



ically and is highly laminar. In this region, which we shall 
refer to as the infall region, we find that the polytropic sim- 
ulation produces quite similar results to the ballistic one, 
and hence, also to the toy model. On the other hand, in 
the post-shock region the flow is decelerated, the stream- 
lines deviate away from the ballistic trajectories due to the 
action of pressure gradients and, in this way, are preven- 
ted from having a 'head-on' collision with their symmetric 
counterparts coming from the opposite hemisphere. Clearly, 
the full hydrodynamical evolution depends on the particu- 
lar EoS being used as well as on the particular mechanism 
driving the accretion (e.g. viscosity, dynamical instabilities, 
etc.) and the cooling prescription, but again, the details of 
this post-shock region were outside the scope of the present 
study. From Figure [4] we also note that the isodensity con- 
tours of the polytropic simulation in the infall region are less 
noisy than those of the ballistic simulation, due to the action 
of pressure forces in smoothing out the particle distribution 
and so reducing discretisation fluctuations. 

Note that, in comparing Figures [5] and [G] the departure 
of the ballistic streamlines from the analytic solution occurs 
earlier and for a larger fraction of the simulation domain 
than the differences between the ballistic and polytropic 
streamlines. In other words we see that here, adopting an 
improved description for the gravitational field of the BH 
has a greater effect on the infall part of the simulation than 
including pressure. 



7.2 Example II 

In order to make a connection with the results presented 
in Paper II and to demonstrate the use of the toy model 
as a useful tool for studying the effect of different gravity 
treatments, we considered the same boundary conditions as 
in Eqs. (|7.2p - (|7.7|) but now with a = 0. This same set of 
parameters was also used in Paper II t o make a comparison 
with one of the GRB simulations of I Lee fc Ramirez- Ruij 
l|200^ using the PW potential. 

Figure [7] shows a comparison of the analytic general re- 
lativistic streamlines with those coming from ballistic SPH 
simulations with the three different gravitational potentials: 
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MM potential 




10 

R/M 
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R/M 



Figure 8. This figure compares the streamlines from SPH sim- 
ulations for ballistic motion against those for a polytropic fluid. 
Note that for just the run using the PW potential the resulting 
flow corresponds to a 'small-scale inviscid disc'. In the other two 
cases, infalling matter keeps accumulating in a ring around the 
BH. This effect is more evident in the case with the MM potential. 



the classical Newtonian one, and the PW and MM pseudo- 
Newtonian ones. Note how the streamlines obtained with 
the Newtonian potential are closer to the general relativity 
solution than those obtained with the pseudo-Newtonian po- 
tentials. In this figure we can also see that the PW stream- 
lines arrive at the equatorial plane at smaller radii than the 
analytic relativistic ones, while the Newtonian and MM ones 
arrive at larger radii. This suggests that equivalent simula- 
tions implementing the PW potential would underestimate 
the extension of any resulting disc while those implementing 
the Newtonian and MM potentials would overestimate it. 

Once more we analyse the role of pressure gradients on 
the infall by showing in Figure [8] the results from polytropic 
simulations for the three effective potentials employed in 
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Figure [7] in comparison with the equivalent ballistic ones. 
We see a quite good match between the two sets of stream- 
lines in the infall region while the effects of the pressure 
gradients become significant only in the high density region 
near to the equatorial plane. Again, in this region we observe 
that a shock front develops surrounding the disc where the 
incoming streamlines decelerate and deviate from the cor- 
responding ballistic trajectories. 

With respect to the a 7^ case in Section rfTI the change 
of the spin parameter of the BH does not lead to signific- 
ant differences in either the velocity or the density fields of 
the accretion flow in the infall region. Nevertheless a cru- 
cial difference comes at the level of the long term evolution 
of the resulting disc since, for the present boundary condi- 
tions with = 0, one expects that none of the matter in the 
equatorial plane would possess enough angular momentum 
to maintain a stable orbit around the central BH and, hence, 
all of the infalling material should be accreted onto the BH 
horizon in a dynamical time-scale. This kind of accretion cor- 
res ponds to the 'small-scale i nvisci d disc' regime discussed 
bv lBeloborodov fc Illarionovl (|200lh . In this respect, an im- 
portant difference among the accretion flows corresponding 
to different effective potentials is already apparent in Fig- 
ure [HI In this flgure we see that only the disc correspond- 
ing to the PW potential evolves as an accretion disc within 
the 'small-scale inviscid disc' regime (in which a stationary 
state is rapidly reached). For the other two cases we observe 
a growing ring of matter with enough angular momentum 
to avoid direct accretion onto the BH (even though the bal- 
listic streamlines for the Newtonian potential appear to give 
the best match to those of the analytic relativistic solution 
in Figure (Tj ; the resulting disc in each of these cases is then 
expected to evolve on a viscous time-scale rather than the 
much shorter dynamical time-scale. 



8 DISCUSSION AND CONCLUSIONS 

In this paper, we have presented an analytic toy model for 
the relativistic accretion of non-interacting particles onto a 
Kerr BH. Taking the assumptions of stationarity, axisym- 
metry and ballistic motion, we have given analytic expres- 
sions for the streamlines and the velocity flelds as well 
as a simple numerical scheme for calculating the density 
field. This model is a generalisation of the one presented 
in Paper II for Schwarzschild spacetime, and it has been 
demonstrated how the earlier results are easily recovered 
from the present solution in the non- rotating limit. 

Using a single analytic expression for describing the 
streamlines constitutes a novel way of expressing the 
solution to the latitudinal and radial motion of timelike 
geodesies in Kerr spacetime. The generality of this expres- 
sion is shown in the Appendix by using standard identities 
of the Jacobi elliptic functions. 

We have explored the effect of frame dragging on the 
resulting accretion flow and found that an effective coupling 
occurs between the BH spin and the angular momentum 
of the infall, leading to more extended discs if the flow is 
co-rotating with the BH and smaller discs in the counter- 
rotating case. 

Our model allows for a fairly wide range of boundary 
conditions to be used, making it an ideal tool for explor- 



ing the effect of different flow parameters (accretion rates, 
angular momentum and density distributions, etc.) in ap- 
plications where the approximations of steady-state and 
axisymmetry are reasonable ones. These assumptions are of- 
ten met in some interesting astrophysical scenarios such as 
under- luminous accretion towards supermassive BHs, wind- 
fed X-ray binaries and coUapsars in which the accretion disc 
remains thin either due to efficient cooling or because it 
evolves within the 'small-scale inviscid' regime. In this pa- 
per we have shown a series of comparisons between the toy 
model and fuU-hydrodynamic, numerical simulations for a 
coUapsar-like setup. Rather good agreement was obtained 
between the simulations and the toy model, under circum- 
stances where one might expect to have agreement. The 
main discrepancies between the resulting accretion flows in 
the infall region have been shown to be related more to the 
different treatments of the gravitational field produced by 
the BH rather than the ballistic description of the infall. 
Indeed, we observed that the effects of pressure gradients 
tend to become important just in an immediate proximity 
of the disc, where a shock front develops and decelerates the 
incoming flow. A new kind of exploratory simulation can 
be envisaged in which simple but general boundary condi- 
tions are set far away from the central object and then, 
by using the toy model, transported down to the region in 
which pressure gradients become dominant where a proper 
hydrodynamical study can then be performed. This kind of 
approach would greatly reduce the spatial domain of the 
simulation, allowing greater resolution and largely reducing 
the computing time. 

Given the analytic nature of the present model, it 
provides a very practical tool for use in benchmarking gen- 
eral relativistic hydrodynamics codes and this, indeed, was 
the main motivation for the present work which forms part 
of a larger project for building a new general relativistic SPH 
code. Also, this toy model allows simple and direct comparis- 
ons of approximate methods for including general relativistic 
effects in simulations on a case-by-case basis. We have used 
it here to test the performance of two pseudo-Newtonian 
potentials (MM and PW) and found that neither is partic- 
ularly well-suited for the off-equatorial motion considered 
here, in terms of either velocity or density profiles. 
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APPENDIX A: 
MOTION 



RADIAL AND LATITUDINAL 



In Section [4] we have given analytic solutions for the radial 
and latitudinal motion of a timelike geodesic in Kerr space- 
time. These expressions are general but involve the use of 
complex quantities in some cases. In this Appendix we con- 
sider these special cases for both types of motion and, by 
mea ns of standard identities for Jacobi elliptic f unctions (see 
e.g. ICavlevlll96l] : lAbramowitz fc StegunI flQTOl ) . we rewrite 
the solutions when necessary in such a way that just real 
quantities are involved. 



Al Radial solution 

Let us consider the solution to the radial integral in Eq. (14. 2p . 
i.e. 



^ 2Vk) ' l + fcsn2(^,fc)- 
Defining u = cn(ip, k) and inverting Eq. ()A.8|) gives 



cn ^(u, k) = — 3_ cjj ^ 
2Vk 



l-fc(l-w^) 1 + k 



l + fc(l-M2)' 2Vfc 



(A.8) 



(A.9) 



from where, and after some algebra, we can rewrite "I>(r) in 
Eq. (1X2)1 as 



$(r) 



1 



I3ra— ara + (a — I3)r 
I3rd + arc — {a-\- 0)r ' 



(A.IO) 



which is now an explicit real function of r. 
Case III: Two pairs of complex conjugates 



dr' 



cE'(r). 



(Al) 



The general solution to Eq. (|A.1|) was given in Eq. (|4.6p . 
That expression involves the use of complex quantities when 
TZ{r) has complex roots and should be handled with care 
when e = 0. In order to give alternative expressions for 
these cases, we consider the various possibilities one by one: 

Case I: Four real roots 

The labelling of the roots here proceeds as described 
in Section 14.11 In this case the solution given in Eq. (|4.6|) 
involves only the use of real quantities. For the sake of com- 
pleteness, we reproduce it here 



This case is a new possibility for Kerr spacetime that is 
not present in the Schwarzschild case since there one of the 
roots of TZ{r) is always zero and, therefore, there is at least 
one other real root (see Paper II). 

We start by noting that, since "l>(r) is defined in 
Eq. HA.1|I in terms of an integral with a complex number 
as its lower limit, "l>(r) is itself a complex function of r. In 
the following, we first split "l>(r) into its real and imaginary 
parts and then show that the latter is independent of r, and 
so <I>(r) — ^{ro) will always be a real function of r. 

According to the labelling of the roots in Section UjT] in 
this case we have r„ = r/ and rt = r/, with Re(r-„) <Re(ri,). 
Also, it is simple to check that a and /?, as defined in 
Eq. (|A.4)I . now form a complex conjugate pair, i.e. a — P* . 
We then introduce the following real constants 



$(r) = 



2cn" 



(rd-ra)(rh-r) 
irt,-ra){rd-r) ' 



iy£(r„ - ra){rd - u) 



(rb — ra){rd — To) 
{rd — rb){ra — r„) ' 



(A.2) 
(A.3) 



C = 



rg + Td 

2 ' 

a + P 



Td 



2i 

a- 13 



2i 



(A.ll) 



Using these definitions, it is simple to check that k^, as 
defined in Eq. (|A.6|l . can also be expressed as 



Case II: Two real roots and a complex conjugate pair 

We have that r„ and r^j (with r„ < \rd\) are the real 
roots while Vt and form a complex conjugate pair. We 
define the following three real constants 



: = Sign(e) \/ (Vd - rt) {rd - , 

13 = \/(r„-rb)(r,-r,), 
p_ {a + pf-ird-r^f 

4:0/3 



(A.4) 
(A.5) 

(A.6) 



From the definition of k^ in Eq. (|A.3p . it is easy to check the 
following relation 



(1-ffc. 
4fc,- 



(A.7) 



Now consider the following identity for Jacobi elliptic 
functions 



^2 + ^2 ' 



(A.12) 



from which it is clear that kr is still a real quantity. On 
the other hand, from Eq. (jA.10|l and the definitions in 
Eq. (|A.11|) . it is simple to check that 



13 r d — a r a + (a — I3)r . C,v — ri{r — /i) 
13 rd + ara — {a + I3)r rj u + ({r ~ jj.) ' 



(A.13) 



i.e. the argument of the function cn ^ in Eq. (|A.10|) is a pure 
imaginary number. Moreover, from Eq. (|A.12[) it follows that 



(A. 14) 



We now consider the following identity for Jacobi el- 
liptic functions 



i(^, k) 



(A.15) 
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where K is the complete eUiptic integral of the first kind. 
Again defining u = cn((/3, k) and solving for ip in Eq. HA.15|) 
gives 



cn ^(u, fc) = i 
k 



1 



(A.16) 

Introducing $(r) as the following real function of r 



$(r) 



(A.17) 



and combining Eqs. (|A.13|I . (|A.14|> . (|A.16|) and (|A.17|I . it is 
simple to check that 



<I>(r) = <l>(r) 



(A.18) 



This last equation explicitly splits $(r) into a real function 
of r and an imaginary constant. From this, it follows that 

$(r) - <l>(r-o) = $(r-) - $(ro) (A.19) 
is a real function of r. 

Case IV: e = 

In this case we have that one of the roots diverges to in- 
finity and so TZ{r) reduces to a third order polynomial. Here 
there are two possibilities for the roots: either the three of 
them are real, or one is real and the other two form a com- 
plex conjugate pair. Appropriate expressions for each case 
are straightforward to obtain from Eqs. HA.2|I and HA.lOp 
after taking the corresponding limit. See Paper II for an 
analogous procedure. 

A2 Polar solution 

We now return to the polar integration discussed in Section 
14.21 and consider the polar solution as given in Eq. (I4.iup . 
That expression is always a real function of 9 but, for some 
values of Q and e, it might involve the use of complex quant- 
ities as intermediate steps. Here we show how to rewrite ^(6*) 
in a way which involves only real quantities. For doing this, 
we consider the following cases: 

Case I: e ^ 

In this case, it follows from Eq. (|4.9|l that Q ^ and so 
all of the quantities involved in Eq. (|4.10|l are real. For the 
sake of completeness, we reproduce the expression here 



cos 9a 



cd" 



cos 9 

cos 9r, 



ke 



(A.20) 



. , „, cos 9a _1 
'^(9) = — — cos 



cos 9 



(A.22) 



\/Q \ cos t 

As noted in Section f42l this result also follows when a = 
0; hence Eq. (|A.22|I is the expression to use in Schwarzschild 
spacetime (see Paper II). 

Case II: e > and Q > 

Here one has that ke is a pure imaginary number. Using 
the following identity for Jacobi elliptic functions 



cd(i/5, k) — cn 



one can rewrite 'if [9) in Eq. (|A.20p as 



cos 9a 



_i / COS y ~ 
: cn I — , fc, 

COS Oa 



with 



kn 



(A.23) 



(A.24) 



(A.25) 



1 - fc2 M Q + ea^ COS* 6»„ 
Case III: e > and Q < 

In this case, we use the identity 

cn((?!), k) = dn (fc fc"^) (A. 26) 

to transform *(&) as written in Eq. HA.24|) into 



/e a cos 9a 



■ dn 



cos 9 1 



cos 9a ke 



(A.27) 



When (5 = then kg — 1 and, since dn{ip, 1) — sech{ip), 
Eq. (|A.27p can be simplified as 



a cos 9, 



■ sech 



cos 9 
cos 9a 



(A.28) 



Case IV: £ = 



where 



In this case, the expressions in Eqs. (|X20)) . (fOU and 
(|A.27|) can be used without further modification. Neverthe- 
less, note that here one has the possibility of reaching the 
polar axis where the polar coordinate is singular. This hap- 
pens when £ — and Q > —so? , and in this case one should 
take 9a = ^ which, although it is not a formal root of the 
equation Q(9) — 0, does constitute a turning point of the po- 
lar motion since here one has that the polar velocity changes 
sign discontinuously every time that the particle crosses the 
polar axis. 



kg — \/—e a? IQ cos 



(A.21) 



Note that when e = then kg — 0, and since cd(</3, 0) 
cos (</?) then Eq. (K20} can be simplified as 
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